#include <iostream>
#include "Meshes/Uniform.hpp"
using namespace Tuna;

typedef TunaArray<double,3>::huge ScalarField3D;

void readFromFile_DX(ScalarField3D& scf1, int number, string filename, int w); 
void writeToFile_Prof(ScalarField3D& scf, int number, string filename, 
		      double dx, int mon, int w); 
void readMatriz(int, string);


double satE[10][63];
double satMat[63][63];

int main() 
{
  double length_x, length_y, length_z, step, start, end; 
  int num_nodes_x, num_nodes_y, num_nodes_z, mon;

  std::ifstream input_cin ("inputTransProf"); 
  input_cin >> length_x >> length_y >> length_z
	    >> num_nodes_x >> num_nodes_y >> num_nodes_z
	    >> start >> end >> step
	    >> mon;
  input_cin.close();

  string modelo, modvis, directorio1, directorio2;
  cout << " Modelo (Lineal o Cuadratica) : "; cin >> modelo;
  cout << " Razon mu_w / mu_o (Mu1 o Mu5) : "; cin >> modvis;
  directorio1 = "./ECLIPSE/sat" + modelo + modvis + "/11saturaciones";
  directorio2 = "./ECLIPSE/sat" + modelo + modvis + "/sat";
  ifstream fileECLIPSE(directorio1.c_str());
  if (!fileECLIPSE) {
    cout << "\n !Can't open file \"" << directorio1 << "\" " << endl;
  } else {
    cout << "\n Opening " << directorio1 << endl;
  }
  int NsatE;
  fileECLIPSE >> NsatE;
  for(int iS = 0; iS < 63; iS++) {
    for(int nS = 0; nS < NsatE; nS++) {      	
      fileECLIPSE >> satE[nS][iS];
    }
  }


  StructuredMesh<Uniform<double, 3> > mesh(length_x, num_nodes_x,
					   length_y, num_nodes_y,
					   length_z, num_nodes_z);
  double dx = mesh.getDelta(X), dy = mesh.getDelta(Y), dz = mesh.getDelta(Z);
  //  ScalarField3D p( mesh.getExtentNodes() );
  //  ScalarField3D S( mesh.getExtentNodes() );
  ScalarField3D p( mesh.getExtentVolumes() );
  ScalarField3D S( mesh.getExtentVolumes() );

  double errorL2, errorL1, errorInf, diferencia;
  int pausa;
  int iteration = start;
  while (iteration <= end) {
    std::cout << "\n ---- Iteration " << iteration ;
    //    readFromFile_DX(p, iteration, "./pres.", 5);
    readFromFile_DX(S, iteration, "./satu.", 5);
    readMatriz(iteration, directorio2);
    
    errorL2 = 0; errorL1 = 0; errorInf = 0;
    for(int i = 0; i < 63; ++i) {
      for(int j = 0; j < 63; ++j) {
	diferencia = satMat[i][j] - S(i+1,j+1,mon);
	if (diferencia > errorInf) errorInf = diferencia;
	errorL1 += fabs(diferencia);
	errorL2 += diferencia * diferencia;
      }
    }
    errorL2 = sqrt(errorL2 / (63 * 63) );
    errorL1 = errorL1 / (63 * 63);
    cout << "\n ERROR L2 (" << iteration << ")= " << errorL2 << endl;
    cout << "\n ERROR L1 (" << iteration << ")= " << errorL1 << endl;
    cout << "\n ERROR INF (" << iteration << ")= " << errorInf << endl;

    //    cout << Sn; cin >> pausa;
    //   writeToFile_Prof(p, iteration, "./pgnu.", dx, mon, 5);
    writeToFile_Prof(S, iteration, "./sgnu.", dx, mon, 5);
    iteration += step;
  }
  
  return 0; 
}

void readMatriz(int number, string filename) {
    ostringstream extnumber;
    //    extnumber.width(w);
    //    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ifstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    for(int i = 0; i < 63; ++i)
      for(int j = 0; j < 63; ++j)
	file >> satMat[i][j];
}

void writeToFile_Prof(ScalarField3D& scf, int number, string filename, double dx, 
		      int mon, int w) 
{   
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }

    int mony = scf.ubound(secondDim) - 1;
    int inicio = 1; 
    int final = scf.ubound(firstDim) - 1;
    double error = 0;
    int numberCol = number / 100;

    cout << "\n Iter = " << number << " Col = " << numberCol << endl;
    for(int i = inicio, iS = 0; i <= final; ++i, ++iS) {
      file << i * dx << "\t" << scf(i, i, mon) << "\t" << satE[numberCol][iS] << endl;
      error +=  (scf(i, i, mon) - satE[numberCol][iS]) * (scf(i, i, mon) - satE[numberCol][iS]);
    }

    error = sqrt(error) / 63;

    cout << "\n Error " << number << "\t" << error << endl;
    file.close();   
}


void readFromFile_DX(ScalarField3D& scf1, int number, string filename, int w) 
{   
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;
    ifstream file (filename.c_str());

    if (!file) {
      cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    char line[256];
/// 
/// by the time, just jump the first 7 lines
///	
    for(int ii = 1; ii <= 9; ++ii)
	file.getline(line,256);
    
    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);

    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
      for(int j = 0; j <= scf1.ubound(secondDim); ++j)
	for(int k = 0; k <= scf1.ubound(thirdDim); ++k)
	  file >> scf1(i,j,k);

    file.close();
}
